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Abstract 



In this paper we consider flat and open universe models containing a mixture 
of cold matter (dust) and radiation interacting only through gravity, with the aim 
of studying their stability with respect to linear scalar perturbations. To this end 
we consider the perturbed universe as a dynamical system, described by coupled 
differential equations for a gauge - invariant perturbation variable and a relevant 
background variable. 

The phase - space analysis of this dynamical system shows that flat dust - radia- 
tion models are unstable, and open models structurally unstable, with respect to 
adiabatic perturbations. For flat models, there are actually three different regimes 
of evolution for the perturbations, depending on their wavelength and the transition 
scale from one regime to the other is determined by a critical wavenumber for the 
perturbations, kEC ( an invariant of the model). We find that k^c < kjE (the 
Jeans wavenumber at equidensity of matter and radiation), implying that there 
are perturbations which decay even if their wavelength at equidensity is larger 
than the corresponding Jeans scale. We also briefly discuss metric and curvature 
perturbations. 

We believe that this analysis gives a clearer idea of the stability properties of 
realistic universe models than the standard one based on the Jeans scale, despite 
our simplifying assumptions. 
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1 Introduction 



A main goal of present - day cosmology is to understand the formation of the structures 
(galaxies, clusters, superclusters) observed in the universe, while trying to explain why, 
on large enough scales, this seems to be so well described by the Friedmann- Lemaitre- 
Robertson- Walker (FLRW hereafter) isotropic models. Given these latter, most theories 
of structure formation are based on the gravitational instability scenario. At any given 
epoch, there are perturbations larger than a certain characteristic - time dependent - 
scaleQ; while perturbations much smaller than this scale oscillate as sound waves, the 
larger density perturbations grow, eventually entering a non- linear regime during the 
matter dominated epoch, thus forming the observed structures. The mathematical basis 
for such scenario is the theory of perturbations of FLRW models. 

In the inflationary scenario, the perturbations were generated from quantum fluc- 
tuations within the Hubble horizon H" 1 , have evolved classically outside the horizon, 
and have re-entered it during the radiation or matter dominated epochs (perturbations 
with larger wavelength re-entering later). During the last decade, observations of the 
distribution of matter have shown that the scale at which the background homogeneity is 
reached is larger than what was thought before, being of the order of hundreds of MpcQ. 
Consequently, both the inflationary scenario and the large scale observations motivate 
the use of a fully relativistic theory of perturbations in FLRW models in order to study 
the formations of the larger structures!] 

For the latter ones, the density contrast against the homogeneous background appears 
to be small enough that a linear perturbation analysis still suffice to describe the evolution 
of perturbations of the corresponding scale. Even if the density contrast is mildly non- 
linear, the curvature perturbations are still in the linear regime ||, [13|], thus one can 



imagine the present universe as well described by a linearly perturbed FLRW model at 
large scales, while non-linearities at smaller scales can be considered smoothed out in 
this picture]^ Then the question arises if in this respect we live in a special epoch - an 
epoch in which large scale perturbations are still in the linear regime - or if this is a 
natural output in the theoretical context of perturbed FLRW models. This question is 
of the same sort as that posed by the "flatness problem" : if the geometry of the universe 



1 This can be the Jeans scale, or the Hubble radius sometimes loosely referred to as horizon; a more 
important scale is actually the sound horizon, see j^j . 

2 The answer to the question "At which scale in the universe there is a transition to homogeneity?" 
depends on how the question is posed, i.e. how such a scale is defined. In a very loose sense we can say 
that this scale is measured by the size of the largest structures we see (of the order of 10 2 Mpc), and 
in a strict sense such a scale does not exist at all if there is not a cut-off in the perturbation spectrum. 
Various reasonable definitions of such a scale can be given (e.g. specifying a level for the fluctuations, 
such that we can say that there is a transition to homogeneity at the scale where the fluctuations fall 
below the specified threshold): see for examples |^6) and |]l0f . 

3 Newtonian perturbation theory being applicable only for vanishing pressure and at scales much 
smaller than the horizon. 

4 In general, this raises the issue of defining a proper averaging procedure, here we do not address 
this problem; see e.g. Q or [fL8| . 
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is non-flat, then we live in a special epoch in which the density parameter is still close 
to unity fi ~ 1. It is well known that in a flat dust model there is a perturbation mode 
that grows unbounded, while in an open dust universe there is a mode that freezes in 
at an epoch z ~ fi _1 . In both cases the perturbation equation ([I]) is the same for all 
modes, because for c 2 s = w = (dust) the coefficient (3 (^3j) does not depend on the 
wavenumber k. It is usually said that in a dust universe each perturbation evolves as 
a separate FLRW universe. On the other hand, perturbation scales in a pure radiation 
model always come within the "sound horizon" (see 0) and oscillate as sound waves. 

This simple picture gives however rather little information about the generality of 
this behaviour for the perturbations, and on the stability properties of the perturbed 
FLRW model. 

The aim of this paper is therefore to study the problem of the stability of FLRW 
models following an alternative approach: instead of looking for analytic solutions (either 
exact or approximate, see e.g. P5f ) of the perturbation equations, we consider the 
perturbed universe as a dynamical system, described by coupled differential equations 
for a gauge - invariant perturbation variable and a relevant background variable. In this 
approach, the evolution of perturbations is represented by the trajectories in the phase 
space of the dynamical system, and their final fate is linked to the presence and the 
nature of the critical points of the system. 

In a certain sense the present work extends that of || to include perturbations, 
although we restrict our analysis to open and flat models only, and we take a vanishing 
cosmo logical constant A = 0. Previous works have followed the approach to stability 



used here p8| , Q, but they only considered either a pure radiation or a pure dust 
model. However, since the analysis here is based on the study of the dynamical system 
at time - infinity, the pure radiation model does not appear physically significant in this 
limit, while the stability properties of the dust models are affected by the simplifying 
assumption of the vanishing of the speed of sound in these models. Therefore here 
we consider a class of perturbed models containing both dust and radiation as a more 
realistic description of the real universe. The dust component can be taken to represent 
Cold Dark Matter (massive weakly interacting particles), while the radiation component 
represents photons and other massless particles such as massless neutrinos. Therefore 
these two fluids interact only through gravity. 

Since our aim is to study the stability properties of these simple models, we consider 
only the total density perturbation (the single component perturbations are not directly 
relevant to the evolution of curvature perturbations), assuming adiabatic perturbations. 

For the case of flat universe models we find that there exists a critical wavenumber 
k^Ci which is an invariant characteristic of the model and is related to the only scale en- 
tering the flat models, i.e. the Hubble radius at equidensity of matter and radiation . 
The corresponding critical scale Xec (-\e is the perturbation wavelength at equidensity), 
remarks the transition from stability to instability, but in a way which is more rigorous 
- from the point of view of the stability analysis - than the Jeans or the Hubble scale. 
We find that there are actually three regimes for the evolution of fluctuations: 
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1. growing large scale perturbations (unstable modes) 

2. overdamped intermediate scale perturbations 

3. damped small scale wave perturbations, 

where the transition scale from one regime to another is always of the order of Xec-, an d 
for this latter, separating cases 1. and 2., we find that the perturbed FLRW models are 
structurally unstable |]. Also, we show that Xec is of the order of the Jeans scale Xje at 
equidensity in the same model: however since Xec — 2-2Xje (today A c — 67Mpc), our 
analysis shows that there are perturbation modes that decay, despite that their scale Xe 
is larger than Xje at equidensity. Thus the evolution of perturbations in these models 
depends on their scale, in such a way that smaller scales evolve like in a pure radiation 
model (case 3.) and larger scales like in a pure dust model (case 1.), while we found 
a small intermediate range of scales (case 2.) for which perturbations are overdamped 
(critical damping occurs for the transition scale between cases 2. and 3.), which is an 
original feature of the dust-radiation models, and to our knowledge was not known 
before. 

For the case of an open dust-radiation model instead, the evolution of perturbations 
appears to be dominated by the curvature of the background, and their final state is 
similar to that they have in a pure dust model, i.e. all the perturbation scales are 
frozen in to a constant value. These models are structurally unstable [[]J with respect to 
perturbations of any wavelength. 

The plan of the paper is as follows: in section ^| we give a brief outline of the method 
used here to study the stability of FLRW models; in section ^| we briefly describe the 
dust - radiation models and we specify the perturbation equations for them; in section 
|4] we describe our results and we compare them with the Jeans stability analysis, also 
we briefly discuss metric and curvature perturbations. Finally, in section ||] we give a 
summary and discuss other possible applications of the method used here. 

In the following we take c = 1, and k = 8irG (G is the gravitational constant), and 
we assume a vanishing cosmological constant A = 0. 



2 The dynamical system 

In a recent series of papers, Woszczyna and colleagues considered the dynamics of New- 
tonian |3(J and relativistic |28]| linearly perturbed universe models, studying the stability 
of these dynamical systems. The analysis of the relativistic case was however affected by 



a wrong assumption on the allowed range for a scale parameter k (q in |28] and in 
and by a misinterpretation of the perturbation variables, as it was shown in H. 

In this article, we shall focus on relativistic perturbations of flat and open FLRW 
models containing two uncouple fluids coupled only through gravity: dust and radiation. 
Open and flat FLRW models expand for an infinite amount of time and therefore one can 
apply standard stability criteria: in particular one can establish if a cosmic dynamical 
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system, i.e. a perturbed cosmological model, is able or not to return to the original 
homogeneous state (in the latter case we have an instability in the sense of Liapunov 
and the system is said to be structurally unstable |L]) or if the perturbations grow 
indefinitely (instability in the Laplace sense). 

In general, in the gauge- invariant approach to cosmological perturbations the scalar^ 
density perturbation variable (we shall consider specific perturbation measures later) 
satisfy a second order (in some time variable) differential equation. If one restricts the 
attention to the harmonic component X of the perturbation, and assumes that this is 



adiabatic (see section pT2|) , its evolution is given by a homogeneous ordinary differential 
equation 

X + a(t)X + p(t)X = , (1) 

where t here is proper time, and the dot indicates a derivative with respect to t. When 
it is not possible to find a simple solution to ([!]), a qualitative analysis of its properties 
is useful in order to determine the late - time behaviour of the perturbations^. The 
coefficients a and (3 are functions given by the background dynamics, but in general 
their time dependence cannot be explicitly determined. Therefore, it is useful to think 
of a and (3 as known functions of one or more parameters, and add to ([!]) the evolution 
equation for the parameters in order to have an autonomous system. A sensible choice 
followed in [[28| is given by the density parameter Q, which in a FLRW model satisfies 
the equation 

n = n(n - i)d + w)e , (2) 

where G = 3H = 3a /a is the expansion of the cosmic fluid, a is the FLRW scale factor, H 
is the Hubble parameter, and w = p/ fi is the ratio of the pressure to the energy density. 
It is useful to change the independent variable from the proper time t to a function of 
the scale factor a; with the choice ofQ r = In a 3 , (^ = 9) equations (jl]) and (H) give 

X" + 4>X' + £X = , (3) 



Q' = Q(Q - + w) 



(4) 



and in general 



w 



:i+w) 



C 2 sW 



W 



(5) 



5 It is standard to call scalar perturbations those related with density perturbations describing the 
clumping of matter (see e.g. |27j]). In this paper we shall consider only these perturbations, as the only 
relevant to the problem of stability of the universe. 

6 One can obviously find numerical solutions, but the study of the phase space allows us to obtain 
general conclusions for a generic set of initial conditions. 

7 There are various typos in psj| : there, the power appearing here in the definition of r is missed in 
equation (3), and consequently a factor 1/3 is missed in equation (9). 
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gives the evolution of w, and in the previous equations the prime refers to the derivative 
with respect to r = In a 3 . 

In a single fluid FLRW model the dynamics is fixed by an equation of state p = wfi, 
with w = const (e.g. w = dust, w = 1/3 radiation), and in this case c 2 = p/fi is the 
speed of sound, and c 2 = w. In dealing with two or more fluids however, w' ^ 0, and 
c 2 7^ w is no longer the speed of sound unless the fluids are coupled, and ip and £ are 
given by 

i/>(w,Q) = ^ + ±-a(w,n) = -\-];(l + 3w)n + ±-a(w,n) , (6) 

Z(w,n) = (7) 

while the second step in (||) is given by the Raychaudhuri equation (k = 811G) 

9 + |e 2 + |/?/i(l + 3u;) = 0, (8) 

governing the evolution of 9 (see e.g. ||). In the simplest case of one single fluid [p8fl , 
[|J with w = const we have a third order autonomous system, given by ([3]) and (ffl) (or 
a corresponding pair of first order equations); also, equation @ forms an autonomous 
first order subsystem in this case. In the most general case we have a fourth order 
autonomous system, with (||) as autonomous subsystem (as we shall see in section j3?T| 
c s = c s( w ) i s fixed once the fluid components are specified). The order of the system can 
however be reduced: first, it turns out that in the practical case (see section |3.1| ) either 
w or Q can be eliminated; second, here we are not really interested in the evolution-law 
of X, but rather in its qualitative behavior. Because of this, we can achieve a further 
dimensional reduction of the phase space passing to the Riccati equation corresponding 
to (pi). Introducing Y = X' and 



U = Y/X = X'/X, (9) 
K = VX 2 + Y 2 , (10) 

we pass in the new phase space {7Z,U, w, Q}, where 

11 u 2 + \ ' 1 ' 

W = -U 2 (12) 

while the evolution of Q and w is still given by (|^) and (|5|) . The variable 1Z represent an 
"amplitude" of the perturbation and it is not directly relevant to the present analysis. 
Since (0) and (|4j), ([5]) form an autonomous subsystem, one can restrict the analysis to 
the phase space {U, w, f2}; moreover, as we said above, in practical cases we can restrict 
our attention either to {U,w} or to {U,Q}. The relevant variable here is U: when it is 



6 



positive we have either a growing density enhancement or an increasing energy deficit, 
while U < indicates that the inhomogeneity is decreasing (note that from its definition 
(|D U is a tangent in the original phase space {X, Y}, thus — oo < U < oo). One is 
therefore interested in the nature of the critical points (if any) on the Q = or the w = 
axis (the final state of the cosmological dynamical system): a stable node on the U > 
semi-axis will indicate that the given perturbation mode will indefinitely grow, thus 
giving a Laplace instability; a stable node on the U < semi-axis will indicate that the 
given perturbation mode will finally decay, i.e. the system is stable with respect to that 
perturbation; a stable node on U = means that the perturbation will asymptotically 
approach a constant value, i.e. the system is Liapunov unstable, but still stable in the 
sense of Laplace (structural instability |]). Finally, if there are no critical points on the 
Q = (or the w = 0) axis, the perturbation maintains a sound wave character at any 
time. 

A discussion of the system ( P^|) (|) can be found in [^J (and reference therein); a 
discussion of the flaws of the application of the analysis given in p8| to the perturbation 
equations of Bardeen and Ellis - Bruni |LT]| , fll2| , (see also f55fl ) is given in j4|. 



3 Dynamics of the dust - radiation models 

We shall now apply the general method outlined in the previous section to the case of 
uncouple dust and radiation. In the following, we shall normalize the scale factor at 
equidensity of dust and radiation, introducing 5* = a/ag. An analysis of the phase- 
space of FLRW models containing dust and radiation has been recently given in (see 
also Jl7j and [p5|| ); here we simply review well known results, that are needed for the 



perturbation analysis, with emphasis on a useful parameterization. 



3.1 The background 

Since the two fluids are uncoupled, we have separate energy conservation, with /i d = 
^(IeS~ 3 and fi r = l/i^S 1 " 4 , where fiE is the total energy density at equidensity. Then 
the total energy density /i = |/is(S' _3 + S~ 4 ) is also conserved, while the total pressure 
is that of the radiation component, p = p r = |/x r = l/i^S 1-4 ; therefore 

W =3(sW (13) 



i j u ' 1. 1 1 i<uju n \ ( 

for dust and radiation we have 

■2 



For two non interacting fluids the quantity (? s = ? is only formally the speed of sound: 

(14) 



s 3(4 + 35) ' 

At equidensity S = 1, w = 1/6, and c 2 s = A/21. Clearly, (TTJ) can be inverted to give 
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and then the expansion of the universe model we consider is parameterized by w, with 
w (1/3 > w > 0) varying from a pure radiation dominated (t — > 0) to a pure matter 
dominated (t — > oo) phase. 

Note that fll3D , with QUI) , is in practice an integral for (|5|); also, given ( |13|) we can 
integrate (|2|): 

n E (s + i) 

n E (s + i) + 2(i-n E )s 2 ' 1 J 

since today So = ao/ a E = 1 + > 1, the density parameter at equidensity is fig ~ 
( 1 + ) • Clearly, equation (|i~6l) can also be inverted to parameterize the expanding 



model with £1. However in the following we shall find convenient the parameterization 
in w, which will allow us to give a unified treatment of flat and open models: from fllSD 
and ([L4D we have 

c ° = 37TT^) ' (17) 



and this with (|5|) give 

w' = w(w-() , (18) 

where hereafter the prime stands for the derivative with respect to r = InS 3 . 
From (fl8|) and (f|) we get 

n = ^ (19) 

which can be used in expressions for £ and ip to obtain these coefficients as functions 
of w only, thus reducing the effective phase space needed for the stability analysis to 
{U,w}, i.e. that of the plane autonomous system fll2|) , dl8|) . We note also that the 



system (|18|), (g), that describes the evolution of the background universe model is a 
plane autonomous subsystem in the full space {U, Q, w}: this will always be a property 
of a dynamical system describing a perturbed universe, since the basic assumption of the 
perturbation analysis is that of neglecting the backreaction of the perturbations on the 
dynamics of the background model. 

Finally, for completeness, we note that also the total energy density is given by 

27 w 3 

" = T^o^W • (20) 

This obviously also follows from the conservation equation // = + to), integrated 

using (|18|); and /i^ooasif^l/3, fi ^ a.s w 0. 



3.2 The perturbed model 

We shall now consider the dynamics of the perturbed dust - radiation models. 

First, we have to specify a perturbation measure X, i.e. the variable appearing in 
([!]) and (|3|). Here, we shall focus on density perturbations, and we shall take A = 
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^ 3 ^V 2 /i)//i as our fundamental variable, as was originally defined in [0 (see also f29| 



ff 



following the covariant approach to perturbation introduced in 

This covariant quantity is an exact measure of density inhomogeneity, it is scalar and 
locally defined. With respect to a FLRW background A is a gauge - invariant variable 
that, once expanded at first order in the perturbations, is proportional to the Bardeen 
variable e m (see ||); thus its components with respect to an orthonormal set of scalar 
harmonic functions are proportional to the density perturbation in the comoving gauge 
(see 0; for a comprehensive treatment of perturbations in this gauge see |2~T| ). 

In the following, we shall restrict our attention to the harmonic components of A, 
where the scalar harmonics Q are defined by 



(3) V 2 Q = -%Q 



(21) 



( 3 ) V 2 is the Laplace operator in the 3-surface of constant curvature and k > for K = 0, 



but k > 1 for K = — 1 [|l(J, ||22|| . Thus in the flat case the wavenumber k is simply related 
to the physical scale of the perturbation i.e. its wavelength A = 2na/k, but this is not the 
case for the open models. However k can always be taken as invariantly characterizing 
the scale of the perturbation, and we shall do so. 

Since the harmonic component A^ of A and e m are just proportional ||, in the 
following the perturbation measure X appearing in (fj) and (§) can be identified either 
with A( fc ) or with e m . 

In the case of a mixture of dust and radiation, in general the evolution equation for 
X is coupled to the evolution equation for an entropy perturbation variable (see fl9|| ; ||; 
[f2~3|1 and also [p5fl). However the coupling is important only at small scales: here we shall 
only focus on perturbations at scales of the order of or larger than the Hubble radius 
H^ 1 at equidensity, thus we shall restrict to purely adiabatic modes, i.e. solutions of ([!]), 
neglecting the coupling with the entropy perturbation. 

In general the coefficients a and (3 in ([[]) are given by (see |jl2fl ; 0; p5fl ) 



a 



(2 + 3c 2 - 6w)H 



(| + Aw - \w 2 - 3c 2 ) k/j, + 12(c 2 - u>) — 



using (|T|]) for the dust - radiation background we get 



+ c — ■ 



(22) 
(23) 



a 



P 



1+w) 



w 2 {5-3wY 



w(l — 3w) K 
(1 + w) ^ 



+ 



Aw k 2 
3(1 + w)a? 



(24) 
(25) 
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From this, we obtain 



i i 1 

iP(w,Q) = -(l + 3 W )(l-0)+ , (26) 

o(l + if J 

6 y 1 + w / 9(1 + io J 

where fc) = 2 7(i+ w ) * s a function that takes a different form, as function of w, 

depending on the curvature K: in an open universe a 2 H 2 = (1 — Q)~ , but in the flat 
case we cannot use the Friedmann equation to substitute for a 2 H 2 . Rather, we use it to 
substitute for 3H 2 = Kfi, and \i is given by (|20|) , while here a = Scle and S is given by 



(Era). Thus 



K = S = %7^S n A; 2 



SlaiHi 



^E E 28 

, 1 _ ( )\ / .- 

27(l+to) ' 

where 3H E — k^e, an d for open models we can substitute for Q = Q(w) from fll^D in all 
the previous expressions. It is clear from (p8|) that in the space {U,w,Q} the function 
H# is not continuous (for k 7^ 0) on the plane Q — 1, except on the line 10 = 1/3, fi = 1; 
this discontinuity of will play an important role in the behaviour of perturbations. 



4 Results 

We shall now summarize the results we obtained from the analysis of the dynamical 
system given by (|T2"|), ( |TB] ) and (^), with coefficients if) and £ given by ( PB"| ) and (|2"7|), first 
restricting to the flat models and then considering the open ones. 

4.1 The flat models 

The flat perturbed models are described by the subsystem (|T2"D, ( |T8"D substituting Q = 1 
in if> (|26"D and £ ( p7|) and using H for X = given in (|28|) . From this latter we see that 
for flat models there is a characteristic scale that, as we shall see, is related to the late 
time behaviour of the perturbations: this is the Hubble radius at equidensity H^ 1 . It is 
therefore convenient to define 

k E = ~\- = 2^ , (29) 

which represents a wavenumber normalized at equidensity: fee = 2rr \e = H^ 1 i.e. 

= 2rr corresponds to a perturbation wavelength that enters the horizon at equivalence 
epoch. 

Stationary points eventually exist on the w — 1/3 and w = axis, with 



Z4 = f(-^±VV£-4^), (30) 
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Point I: Unstable Node 
w = 1/3, W_ = -i 


Point II: Saddle 

w = 1/3, U+ = | 


X u = l 

\ _ l 

^u) 3 


A w = -1 

X - l 

^to 3 


Point III: Saddle 
w = 0, W_ = - X u 


Point IV: Stable Node 
u> = 0, W+ = - X u 


\ _ l 

3 


x _ l 

Aw ~ 3 



Table 1: Flat models: the four stationary points with the corresponding eigenvalues 
along the axis U and w and their nature; Point III and IV exist only for Ue < itj^ec- 



where ip w , £ w (with w — 0, 1/3) correspond to the stationary values of ip and £. 

The stationary points for the system (|T2|), ([T8|) are given in Table |1|. 

As we have explained at the end of section |2|, the condition for having growing 
perturbations is given by the appearance of a stable node on the positive side of the U 
axis. From (|30D we see that we need £o < in order to have U + > 0, i.e. 

k E < k E c , k EC = • (31) 

Moreover, stationary points only exist for ijj^ — 4£ w > 0, which is always satisfied only 
for w = 1/3. For w = we have — = |(|| — k E /k EC ), i.e. stationary points 

exist on this axis for perturbations of wavelength at equivalence A^ > ^^-Xec < Xec- 
We see from Table |I] that this is the same condition for the reality of the eigenvalues, a 
result that follows from the fact that w' does not depend on IA. From the values of these 
eigenvalues we have that: Point I is an unstable node, Point II is a saddle, Point III is a 
saddle, Point IV is a stable node. 

When for this latter we have U + > 0, i.e. for kE smaller than the critical wavenumber 
kEc, generic perturbations with physical wavelength A# > Xec grow unbounded, with the 
critical perturbation scale Xec corresponding to k E c m (|3ll) , given by Xec — 'i^^E 1 , 
i.e. Xec > H% : this situation is illustrated in Fig. |l|. Note in this figure the special 
saddle trajectories: the one ending in U- represents the purely decaying mode, while 
the other (ending in U + as the generic trajectory) represents the purely growing mode 
(the attractor); the generic trajectory represents a linear combination of the two modes 
solutions of (0). Also, we remind that the variable U is a tangent in the original phase 
space {X, X'}, so that the trajectories in the figures that start from Point I and go to 
the left of the saddle trajectory exit the figure on the left boundary, and re-enter it from 
the right ending in Point IV. 
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Figure 1: Flat models: phase space for the evolution of large scales \e > Xec that grow 
unbounded (U+ > 0). 



Figure 2: Flat models: phase space for the evolution of the damped waves on small scales 

(Z4e9). 
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Figure 3: Flat models: phase space for the over damped modes (U + < 0) on intermediate 
scales ^-Xec < X E < X E c- 

Clearly, when Point IV does not exist, i.e. for perturbations with wavelength at 
equivalence Xe < "^t^Xec-, we are considering perturbations that always oscillate as 
"sound waves" (see Fig. 0); we have checked through a direct stability analysis of the 
system ( |37D below (see section [4.4j) that the amplitude of these modes decay. Again, 



trajectories that exit from the left re-enter from the right: for Xe < ^0-Xec however this 
happens many times, and indeed it is this that characterizes the oscillatory behaviour of 
these modes in the {U, w} plane. 

Also, for scales Icec < &e < j^/q^ec we have U + < for Point IV, i.e. the stable 
node is located on the negative side of the U axis: thus perturbation scales in this range 
are overdamped:they decay without oscillating, as it appears from Fig. Again, the 
saddle trajectories represent the two modes for (|3]), and the trajectories going out from 
the left of Fig. ^| re-enter from the right, ending in Point IV. 

Finally, we point out that for kE = ^^^ec the two stationary points on the U axis 
(the saddle and the stable node in the second line of Table [I]) coincide, i.e. this is a fold 
bifurcation point (see e.g. |1[|). The corresponding scale Xe = ^^-Xec = ^^^H^ 1 ~ 
4.7 H^ 1 is quite larger than the Hubble radius at equidensity; thus we can consider 
scales Xe < ^^-Xeci as m Fig- PL such that still A# > H^ 1 : this partially justifies our 
assumption of purely adiabatic perturbations (cf. 

Thus we have three different evolution regimes for the adiabatic perturbation modes 
of a mixture of uncouple dust and radiation in a flat universe, depending on their wave- 
length: 

1. large scale perturbations that grow unbounded, giving instability 

2. intermediate scale perturbations that are overdamped, i.e. decaying without oscil- 
lating 



13 



Instability 


Stability 


^ > ^EC 

grow unbounded 


^^■^EC < ^E < ^EC 

critical damping 


*E < ~^-^EC 

damped waves 



Table 2: Flat models: summary of the three different evolution regimes for different 
perturbation wavelengths. 



3. small scale damped perturbations which oscillate like sound waves while their am- 
plitude decays, 

as summarized in Table 0. 
4.2 The open models 

The open models are in principle described by the full 3-dimensional system given by 
(|T2"D, flT8| ) and (TJ), with trajectories in the {U, w, Q} space. However, given Q = Q(w,Qe) 



(pj|), we can substitute for Q in the expressions for ip (p6|), £ (j27Q and 5 for K — — 1 
(|28|) , and restrict our analysis to the 2-dimensional system (|i~2"l), ([18|) . In doing this, 
we are selecting one particular open model in the class parameterized by Qe'- there is 
no loss of generality in doing this, as the dynamical properties of the models in this 
class (as specified by the stationary points in their phase space) do not depend on O^, 
i.e. for a given wavenumber k the phase space evolution of all the models in the class 
is qualitatively the same. From the point of view of the geometry of the phase space 
{U, w, Q} we are looking at the trajectories in the 2-dimensional surface specified by Qe 
in this space: we shall then consider the projection of this surface with its trajectories 
in the {U, w} plane, as depicted in Fig. £|. 

It is clear from this picture and Table ||] that, contrary to what we have seen for the 
flat models, the existence, position and the nature of Points I-IV do not depend on the 
wavenumber k, so that all open models share the same dynamical history. This depends 
on the vanishing of the function H in the limit w — > 0, which also implies Q — > (because 
we are moving on the surface specified by Qe)- Note however that, if we do not substitute 
for Q from ( |T9| ) in the 3-dimensional system, and we take the limit Q — > 0, this does 
not give the 2-dimensional system for flat models: this fact is due to the discontinuity 
(remarked at the end of section ^4.1| ) of the function S on the surface Q = 1 , where S for 
K = do not vanish in general, except for w = 1/3 or k = 0. In some way, this fact can 
be seen as an example of " fragility" in cosmology . 

From table and Fig. || we see that Point IV is a stable node located at U = 0. The 
generic trajectory ends up in this point, either directly from Point I, or first going out 
from the left boundary of the picture and then re-entering from the right. The saddle 
trajectory starting form Point I and ending in Point III represents the evolution of the 
purely decaying mode, and the saddle trajectory starting from Point III and ending in 
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Figure 4: Open models: the phase space evolution is independent of the wavelength, and 
all perturbations tend to a constant value (U + = 0). 

Point IV represents the purely growing mode^] (attractor). The fact that Point IV is 
located in U = means that all the perturbation modes evolve up to a constant value, 
and then freeze-in. We can say that open models are Liapunov unstable (the generic 
perturbation modes do not decay), but Laplace stable (the perturbations do not grow 
unbounded), a situation that can be described as structural instability of the dynamical 
system fl|]. 

A direct analysis of the system (pTf ) below (section [O]) shows indeed that the appear- 
ance of a fixed point on U = for the system ([12]), ( |i~8D corresponds to the vanishing of one 
of the eigenvalue for the corresponding fixed point in the phase space {X, X', w}. While 
this happens only for fc^ = fc^c m flat models, it is a generic characteristic for any k in 
open universes. A comparison of Fig. [| and Fig. la in |4[] shows that open dust-radiation 
models and open pure dust models share the same stability properties. It appears that 
the curvature of the background dominates also the evolution of the perturbations: for 
dust radiation models there is a conspiracy between curvature and pressure, such that 
curvature has an opposite effect on small and large scale: with respect to the behaviour 
in flat models, it avoids the oscillation of small scales perturbations (Ag < ^^-Xec) and 

damps the growth of those on large scales(A£ > 



8 The terminology growing and decaying is purely conventional: it is adopted here because it is 
standard in the literature to refer in this way to the two modes, as they are effectively growing and 
decaying e.g. in flat pure dust models. We stress again that the "growing" mode is actually growing 
only when 14+ > 0. 
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Point I: Unstable Node 
w = I U- = -\ 


Point II: Saddle 

ty = I, U + = \ 


A w = l 

\ _ l 

"u> 3 


A w = -1 

\ __ l 

"to 3 


Point III: Saddle 
w = 0, U- = -\ 


Point IV: Stable Node 
w = 0, U+ = 


Xu—\ 
X - i 

"to 3 


Aw = — \ 
\ _ i 

"li> 3 



Table 3: Open models: the four stationary points with the corresponding eigenvalues 
along the axis U and w and their nature. 



4.3 Comparison with Jeans instability 

We can apply the Jeans instability criterion directly to equation (0): this implies that 
gravitational collapse of a given perturbation mode will occur if j3 < 0, i.e. if k < kj, 
where for the dust-radiation models j3 is given by fl25|), and the corresponding kj is 



*3 



3(1 +w) 
Aw 




(5-3w) \ + Ml - 3w) K 



1 + w 



(32) 



At equidensity (w = 1/6) this gives 



a% 



279 
^2~ 



0, 



2 . 

E i 



(33) 



in the real universe Qe ~ 1 and the contribution to /cje from the curvature term (the 
last in the equation above) is completely negligible. Thus in a flat universe (K = 0) we 
have 

2naE 87r / 2 
~ ~3~ 



Aj_b = 




(34) 



/cje 3 V 31" 

for the Jeans scale at equidensity. Then a comparison with the critical scale Xec defined 
gives 



in section f4~ 



A 



.IE 




CE 



(35) 



i.e. for flat models the stability criteria used in section |4.1| give a critical scale Xec f° r 
instability which is larger than the corresponding Jeans scale at equidensity by a factor 
of 2, i.e. Xec ~ 2.3 Xje- The fact that the values of these two scales are relatively close 
appears physically meaningful, and in a way obvious, since both scales are somehow 
defined through the same differential equation. However, a comparison of (|3~5|) with the 



analysis of section [O] shows that there are perturbations with A^ > Xje that decay: 
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those with wavelength Xje < Xe < ^^-Xec are damped oscillation, while those with 
Xje > Xe > are overdamped. It is immediate to show that today the critical 

scale corresponding to k E c is 



X c = - ] J-(1 + ze)-'H \ (36) 
i.e. about 67.3/i _1 Mpc in a flat universe with I + ze ~ 4.310~ 5 /i~ 2 (see e.g. jSU|). Even 



if this value may be an artifact of our simplifying assumptions, e.g. the fact that we have 
neglected isocurvature modes at all times and we have treated radiation as a perfect 
fluid even at small scales, we believe that the discrepancy between our critical scale and 
the Jeans scale in the same universe model is a general feature that deserves further 
investigation in order to consider possible effects for models of structure formation in the 
universe. 

Finally, it is interesting to consider the limit w — > of kj. For flat models, one 
gets from (^) that in this limit kj = ^po^if^, i.e. we recover k E c (|31|) in this limit. 
However the same limit for open models gives a value for kj that: a) is real only for 
Qe > 16/25, a result that, although satisfied in the real universe, appears spurious for 
the theory; b) gives the false impression that there could be growing and oscillating 
modes also for open models, contrary to what we have shown in the previous section. 

4.4 Metric and curvature perturbations 

In the previous sections we have given the results of the analysis of the dynamical system 
{U,w} (f[|), (0) for flat and open models, and inferred conclusions on the evolution of 
density perturbations, represented by the harmonic component X (see section \5.'2\). As 



we have already pointed out, there is a particular relation between the location and 
character of the critical points in the phase space {U,w}, and the character of the 
corresponding point in the original phase space {X, X',w}^. Indeed, it is immediate 
that the original system (Y = X') 

X' = Y 

Y > = -ipY-ZX (37) 
w' = w (w — IJ 

admits only two fixed points for £ ^ 0: Point A = {X = 0,Y = 0,w = 1/3} and 
Point B = {X = 0, Y = 0, w = 0}. Then it is easy to see that the eigenvalues at these 
points along the principal directions in the {X, Y} planes equate the roots U± of (|1"2"D, 
i.e. A± = U±. Point A is the same for flat and open models (A_ = —1/3, A + = 2/3, 
X w = 1/3), while (27), (p5[) and the analysis in section |4.2j shows that for open models 



£ = for w = 0, and thus Point B degenerates into a line (the Y = axis): A + = U + = 

9 Here we assume that, for open models, we are using the function £1 = Q(w, Q) (|l9[), so that the 
further dimension in the phase space is suppressed. 
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Quantity 








Xe — A^c* 


Xe 3> Aec 


x ± 


S -i±*Q 


s-i 




S~? , const. 




<$>N± 




5-4 


S~i ±p 


S~-2, S- 1 


iS 1 " 5 , const. 


E ± /e 2 






s -i±p 


, const. 


S--2 , S 



Table 4: Flat models: asymptotic behaviour for w — > (i.e. around Points III— IV) of 
X, $jv and E/Q 2 , in order of increasing wavelengths, as function of the scale factor S 
(S ~ ts for w — > 0). The decaying (— , Point III) and growing (+, Point IV) modes 
are distinguished either by the ± or presented in order: they coincide for the critical 



damping scale Xe 



2V6 



ec; Q=[l (k 2 E /klc ~ 



and P 



(if k E /k 



2 

EC 



in this case, and we have structural instability M. For flat models, the same happens for 
k-E = k E c-, as already pointed out. 

Having clarified the relation between the roots U± of the system dl2|) , (18) and the 
eigenvalues at the fixed points of (0), we can now turn to the asymptotic evolution of 
X. It is clear from the definition of U = Y/X that around the roots U± the evolution of 
X is given by 

X' = U±X , => X± ~ S 3 "* . (38) 

When U + = we have X + = const, (cf. Table |^); for w->0we recover the well known 
constant mode for matter dominated open models: here the same mode is found for the 
critical scale k E = k E c in flat models. 

In tables |4] and [5] we give the asymptotic behaviour for X, for a metric perturbation 
$at and a dimensionless curvature perturbation scalar E/Q 2 . In particular, in Table f| 
we consider flat models in the limit w — >• 0, giving the asymptotic behaviour of various 
variables as functions of the scale factor S and in order of increasing wavelength Xe- 
In Table |5] we give the asymptotic solutions of open models around Points LTV: as for 
w — > 1/3 the universe is radiation dominated and also ft — > 1, then S ~ tz in this limit, 
and the asymptotic solutions around Points I and II are in common with flat models. 
In the limit w —>■ the universe models are matter dominated, and then in flat models 

2 

S ~ £3, while in open models S ~ t. 

In the following we shall outline the relation between the gauge-invariant metric 
potential $at, the curvature variable E, and the density perturbation A (see section 
|3.2| ): more details can be found in || and references therein . Let E a ^ be the electric 
part of the Weyl tensoiQ C ac bd- E ab = C ac bd,u c u d , where u a is the 4-velocity of matter; 
then to first order 

a m v m V a Eab = ^ A ; (39) 
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Latin indices are 4-dimcnsional (0, 1, 2, 3), and greek indices 3-dimensional (1, 2, 3). 
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Quantity 


Point I 


Point II 


Point III 


Point IV 


X 


s- 1 


S 2 




const. 


<$> N 


s- 3 


const. 


s- 2 


s- 1 


E 
02 


s- 1 


S 2 


s- 2 


s- 1 



Table 5: Open models: asymptotic behaviour around Points I-IV for X, $at and E/Q 2 , 
as function of the scale factor S. For w — > 1/3 S ~ tz, and Point I and II are in common 
with flat models, while for w — > S ~ t. Points I and III correspond to the decaying 
mode, and II and IV to the growing mode. 

where ^ V a is a covariant derivative orthogonal to u a . The analogue of E a b in Newtonian 
theory is the tidal field E a p = V Q/ 30, where is the Newtonian potential and V a p(p = 
V 2 — !<5 a/ 30. Then it is possible to show that to linear order a similar formula holds in 
relativistic perturbation theory, i.e. E a p = Va/s^iv for the scalar part of E a b (all the {0, 0} 
and {0,a} components are second order). Moreover the field $tv, which play here the 
role of a gauge-invariant analogue of the Newtonian potential, is just $n = §($a — &h), 
where $^ and are the gauge invariant metric perturbations defined in ]2[] ($^4 = — 
for perfect fluids). Then using E a/3 = V a/ 3^N in (|39|) , the harmonic decomposition (|2l| ) 
and = —k 2 e m , we get equation (4.3) of 

2(3K - k 2 )$ N = Ka 2 fie m . (40) 

Then this relation can be used to determine the asymptotic evolution of the gauge- 
invariant metric potential $tv, and also that of amplitude E = ^y/E a b Eb a of the tidal 
field Eat,: indeed from a comparison of ( |3"9"D and (^) the harmonic components of E and 
$tv are related by 

E = \a- 2 k 2 $ N . (41) 

Then it is usual to consider a dimensionless scalar to measure the relative dynamical 
significance of a given field using to take into account the expansion (see e.g. flT5l ): 
thus in the case of E we consider E/Q 2 , for which we have 

E a~ 2 $N kll „ , 

e5 ~ ~ ^ ~ ■ ( 42 ) 

Hence, for example for a flat model $tv ~ const, for very large scales, but the relative 
amplitude of the tidal field E/Q 2 ~ A ~ a grows unbounded (cf. @). 



5 Conclusions 

In this paper we have considered the stability properties of FLRW models with uncouple 
Cold Matter (dust) and radiation. We have considered the perturbed models as dynami- 
cal systems described by an evolution equation for a gauge-invariant density perturbation 
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variable coupled with the equations governing the evolution of relevant background vari- 
ables: the pressure-energy density ratio w = p/fi and the density parameter Q. For the 
subsystem describing flat models, given by O = 1, we deal with a planar autonomous 
system, and for the open models we have shown that we can also restrict the analysis to 
a planar system for each particular value of the density parameter at equidensity Jig. 

The analysis of flat and open models gives different results: flat models admit unstable 
perturbation modes, while open models are Laplace stable (but Liapunov unstable) with 
respect to perturbations, irrespective of their size, ie. the perturbations freeze-in at a 
constant value (structural instability). Qualitatively, both results can be expected on 
the basis of the results of a more traditional study of the behaviour of the perturbations 
in simple models. The final fate of the perturbations in open models appears to be 
dominated by the background curvature, which governs the background expansion at 
late times: therefore the stability properties of the dust-radiation models are the same 
as those of a pure dust model, irrespective of the size of the perturbation and of the 
radiation content of the given model. 

Instead, we find more interesting features for flat models: in this case each model has 
a characteristic critical invariant wavenumber kEc which depends on the proportion of 
matter and radiation (e.g. at present), and the corresponding scale Xec determines the 
transition from stable to unstable modes. The scale A ec is of the order of the Hubble 
radius H^ 1 at the equidensity of matter and radiation, H^ 1 being the only scale entering 
the background model. The present value of this transition scale is Aoc = (1 + Ze)Xec 
and depends on the actual present proportion of matter and radiation, for which it is of 
the order of 67Mpc. Therefore the stability properties of dust-radiation flat models are 
a mixture of the properties of pure dust and pure radiation models, being well known 
that perturbations in a flat dust models grow unbounded irrespective of their size, while 
all perturbation scales in a radiation model always enter an oscillatory (sound waves) 
regime after they enter the sound horizon (see e.g. @). 

However, we actually find a structure in the stability properties of flat dust-radiation 
models which is more interesting than expected, because we find three different regimes 
for the evolution of the perturbations. Perturbations on scales Xe > Xec grow un- 
bounded (unstable modes), while perturbations in the range ,. A ec ^ A^; < Xec are 

overdamped; finally, perturbations on scales Xe < ^^-Xec are damped, i.e. they oscillate 
as "sound waves" while decaying. 

The stability properties of perturbed FLRW models that we have found should be 
taken cum grano salis with respect to the problem of structure formation in the universe, 
because we have assumed adiabatic perturbations at all times, neglecting isocurvature 
modes that are in principle important at small scales and late times, and we have ne- 
glected photon diffusion, treating radiation as a perfect fluid irrespective of the scale of 
the perturbation (in particular, the oscillatory behaviour of small scales in flat models 
is probably an artifact of this assumption). However, we believe that these assumptions 
should not question the validity of our main results: a) in a given universe model (includ- 
ing the assumptions about the matter content and how to treat it) there are perturbations 
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that decay despite the fact that their wavelength at equidensity is larger than the 
Jeans scale Xje] b) the today size of the critical scale we have found is of the order of 
67 Mpc, a fact that perhaps deserves further investigation regarding its implication for 
structure formation in the universe. 

Finally, we remark that a comparison of the stability properties of dust -radiation 
models with the observed small amplitude of the large scale density fluctuations seems 
to suggest that if the spatial curvature of the universe vanishes, then we live in a special 
epoch in which these perturbations are still in a linear regime of growth, while if we live in 
a universe with negative spatial curvature the smallness of the large scale perturbations 
is a characteristic of the model at all times and the description of the universe as an 
open FLRW model is appropriate at any epoch. 

Thus open models are special (at a given time) on average, because the density 
parameter Q depends on time (this is the flatness problem), but not from the local point 
of view, because large scale structures are frozen-in, while the reverse is true for flat 
models, because Q = 1 at all times, but large scale perturbations grow unbounded in 
these models. 
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